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ABSTRACT 

I argue that the problem of electromagnetically driven electron-positron cascades in mag- 
netospheres of neutron stars must be addressed starting from first principles. I describe a 
general numerical algorithm for doing self-consistent kinetic simulations of electron-positron 
cascades - wherein particle acceleration, pair creation and screening of the electric field are 
calculated simultaneously - and apply it to model the Ruderman & Sutherland (1975) cas- 
cade in one dimension. I find that pair creation is quite regular and quasi-periodic. In each 
cycle a blob of ultrarelativistic electron-positron plasma is generated, it propagates into the 
magnetosphere leaving a tail of less relativistic plasma behind, and the next discharge occurs 
when this mildly relativistic plasma leaves the polar cap. A short burst of pair formation is 
followed by a longer quiet phase when accelerating electric field is screened and no pairs are 
produced. Some of freshly injected electron-positrons pairs get trapped in plasma oscillations 
creating a population of low energy particles. The cascade easily adjusts to the current density 
required by the pulsar magnetosphere by reversing some of the low energy particles. Each 
discharge generates a strong coherent superluminal electrostatic wave, what may be relevant 
for the problem of pulsar radioemission. 
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1 INTRODUCTION 

Rotation powered pulsars remain a profound puzzl e despite the 
fact th at the first pulsar was discovered 40 years ago (iHewish et al.! 
[l968|). Pulsar is a rapidly rotating strongly magnetized neutron star 

i NS), as it was originally proposed by Gold (1969) and Pacini 
1963), with most of its radiation produced in the magnetosphere. 
However, there is still no consistent quantitative pulsar model. Pro- 
posed models range from a NS w ith charge starved electrosphere 
(iKrause-Polst orff & Michef 1985') to a NS with force-free magne- 
tosphere, where acceleration of particles and, h ence, emitting zones 
are lo calized in very small spatial regions (iGoldreich & JulianI 
[T969b . 

The force-free magnetosphere model is favored by majority of 
astrophysicists working on pulsars. There are observational hints 
favoring this model: i) young pulsars produce relativistic winds 
with particle number density much larger than it is necessary to 
screen accelerating electric field parallel to the magnetic field; ii) 
pulse peaks are narrow, what points to smallness of emitting re- 
gions, and, hence, to smallness of regions where particles are ac- 
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celerated . From theoretical point of view, as it was pointed out by 
ISturrockl(ll97lh . physical conditions in the polar cap of pulsar are 
almost ideal for generation of electron-positron pair plasma. The 
energy density of the generated plasma is negligible compared to 
the energy density of the magnetic field near the NS. The magneto- 
sphere, if filled with plasma, almost certainly being force-free (al- 
most everywhere) near the NS should be force-free at much larger 
scales as well; at least numerical simulations of the force-free mag- 
netosphere of an aligned rotator have shown that magnetosphere 
can remain for ce-free up to dis tances much larger than the light 
cylinder radius (iTimokhinlbOOi) . 

Therefore, pursuing the force-free model as a 'standard 
model' seems to be reasonable. Recently the force-free pul- 
sar magnetosphere model has been studied in great detail 



fe.g. IContopoulo s et al.l '19991; 'Oruzinov l l2005l: iTimokhii 
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IBai & Spitkovskvll201 0). Force-free magnetosphere is a restricted 
MHD system which does not admit any current density distribu- 
tion. By fixing the boundary conditions - in the case of pulsar these 
are the variation of accelerating potential across the polar cap and 
the size of the corotating zone - one fixes the current density dis- 
tribution. It turns out that the admitted range of current density dis- 
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tributions in the force-free magnetosphere with reaUstic boundary 
conditions - when the potential drop i n the polar cap is le ss than 
a vacuum one - is rather limited ( TimokhinI l2006l . l2007bl fal). For 
young pulsars, where potential drop in the polar cap must be small, 
the current density is not constant and strongly deviates from the 
Goldreich- Julian (GJ) current density j^^ = rj^jC (t/gj is the GJ 
charge density, c is the speed of light); along some magnetic field 
lines it has the sign opposite to the sign of the GJ charge density. 

Pair production in the polar cap of pulsar is vital for sus- 
taining of the force-free magnetosphere - without it there will be 
not enough plasma to cancel the accelerating electric field. Cur- 
rents flowing in the open field line zone of the magnetosphere flow 
through the pair-producing region at the base of the polar cap; 
therefore, any model of the polar cap cascade zone must agree 
with a global magnetosphere model on the current density flow- 
ing along magnetic field lines. Many previously proposed mod- 
els for polar cap cascades (and almost all quantitative models) as- 
sumed station ary unidirectional outflow of a charge separated parti- 
cle be a m (e.g.lArons & Scharlemann 19791: Daughertv & Harding 



19821: iMuslimov & TsvganI Il992l : faarding & M uslimov I2OO2L 
Hib schman & AronsI 12001 ah . All these models predicted current 
density being almost equal to the GJ current density everywhere 
in the polar cap of pulsar. This prediction is in strong disagree- 
ment with the force-free magnetosphere model: for young pulsars 
like Crab a deviation of the charge density from t/gj of the order 
of few per cents - and in unidirectional flow this implies the same 
deviation of the current density from Jgj - can account for all pul- 
sar emission. Both sides of this discrepancy are based on detailed 
simulations and it is not possible to change some parameters in or- 
der to fit the models together. So, either the magnetosphere is non- 
force-free or non- stationary (or both) or polar cap cascades do not 
operate according to the existing models. 

From the energetic point of view a stationary (on the rotation 
time scale) force-free configuration seems to be the most preferable 
state of the magnetosphere. The inductance of the magnetosphere 
is much larger than that of the polar cap, therefore, the current den- 
sity in the polar ca p will be set by the magnetosphere and not in the 
opposite way (e.g. lMestelll999b . In my opinion these are strong 
hints that existing quantitative models for particle acceleration and 
pair production in pulsar polar cap do not work. Particle accelera- 
tion and electron-positron pair production in cascade zones can be 
essentially non- stationary: time intervals of eff'ective particle accel- 
eration could alternate with intervals when the accelerating electric 
field is screened by electron-positron pairs created in the cascade; 
in fact, in the first paper on pulsar cascades (Sturrock 1971) the 
particle flow was assumed to be non- stationary. The current den- 
sity flowing through non- stationary cascade fluctuates strongly and 
the amplitude of the fluctuation should depend on the microphysics 
of the pair-generation process, not on the global physics of the 
magnetosphere. However, the characteristic timescale of polar cap 
cascades (microseconds) is much shorter than the magnetospheric 
timescale (longer than milliseconds) so that all fluctuations due to 
cascade non-stationarity will be washed out. The average current 
density in the cascade zone could be adjusted to the current density 
required by the magnetosphere by adjusting the time cascade spent 
in "active" and "passive" phases. On the other hand, it is still possi- 
ble that particle flow in the cascade zone is nevertheless stationary 
but not unidirectional - with some particles trapped in a non-trivial 
accelerating potential (Aronsi l2009h . However, all these qualitative 
statements have to be proved. 

Electromagnetically driven electron-positron cascades can op- 
erate not only in polar caps of radiopulsars. Some pulsars in 



outer parts of their magnetospheres - close to the place where 
the GJ charge density changes the sign - co uld have so-called 
"outer-gap" cascade zones (Cheng et al. 1976); although it seems 
that such acceleration zone can exist only if polar cap cascades 
fail to supply enough pair plasma to short-out the electric field 
in the entire magnetosphere. Electromagnetically driven cascades 
should generate plasma in magnetospheres of magnetars along 
open (Thompson 2008) as well as closed magnetic field lines 
(Beloborodov & Thompson 2007). Electron-p ositron cascades ca n 
also work in magnetospheres of black holes jBeskin et al.lll992h . 
The study of pair cascade dynamics is, therefore, of significance 
for a broad class of astrophysical problems. 

Non- stationary regime of electromagnetically driven cascades 
is poorly investigated. Only few attempts have been made be- 
fore to const r uct qu antitative models for non- stationary cascades. 
lAl'Ber et"aD (Il975h were the first, their model was OD - it ac- 
counted only for variability in time. It predicted strong time 
variability in pair creation rate due to the delay between emis- 
sion of a hig h energ y phot on and its decay into an electron- 
pos itron pair. iFawlevI (Il978h tried to make a numerical model 
for lRuderman & Sutherland! (Il975') cascade using ID Particle-In- 
Cell (PIC) code and a simple version of on-the-spot approxima- 
tion for pair injection. At that time it was a formidable numer- 
ical problem; simulations could be performed only for a very 
short time after cascade igni tion, so that no c o nclusive results 
could be dr awn from them. [ Levinson et al.l (l2005h : lLuo & Melrosd 
(2008); M elrose et all (I2OO9I) used ID two-fluid approximation for 
electron-positron plasma and on-the-spot approximation for pair 
injection; they studied polar cap cascades operating in the space 
charge limited flow regime and found that generation of pairs is 
essentially turbulent - pair we re created throughout all physica l 
region admitting pair creation. iBeloborodov & Thompson! (l2007h 
studied pair cascades in the closed field line zone of magnetar mag- 
netosphere. They used on-the-spot approximation for pair injection 
and tracked motion of electrons and positrons in self-consistently 
calculated electric field; the electric field was assumed to be zero at 
both ends of the field line. They too concluded that pair creation is 
turbulent. 

In all of these models some or other simplifying assumptions 
about physical processes at play were used. It is difficult to draw de- 
cisive conclusions about the character of particle flow pattern from 
them because it is not clear a priori whether ignoring one of the 
aspect of cascade physics can result in qualitatively difl'erent be- 
havior or not. In my view, the study of electron-positron cascades 
should be done starting ab initio. No assumptions about the char- 
acter of particle flow should be made and the key "ingredients" of 
the system must be preserved in the model: back reaction of parti- 
cles on the accelerating electric field and the delay between photon 
emission and pair injection. Possible complexity of system behav- 
ior compels to conduct a numerical experiment where particle ac- 
celeration, pair production and variation in the accelerating electric 
field are modeled self-consistently. 

With this paper I intend to start a series of publications ded- 
icated to self-consistent numerical modeling of full kinetics of 
electron-positron pair cascades in magnetospheres of neutron stars. 
In this paper I describe a numerical algorithm for self-consistent 
modeling of electromagnetic cascades starting from first princi- 
ples and apply it for study of the most simple model of polar 
cap cascade - when part i cles c annot escape the NS surface - the 
iRuderman & Sutherland! (Il975h model. The goal of this work is 
not merely to quantify the Ruderman-Sutherland model but to try 
to infer basic properties of electromagnetic cascades. The most im- 
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portant qualitative questions about basic cascade properties I try to 
answer are i) what is the character of plasma flow and ii) how the 
pair cascade adjusts to the current density required by the magne- 
tosphere. 

The structure of the paper is as follows. In Sec. [2] I describe 
the general numerical algorithm I developed for modeling of elec- 
tromagnetic cascades. In Sec. [3] I describe physical and numerical 
aspects of the of polar cap cascade model. Simulations results and 
their analysis are presented in Sec. ID I summarize the inferred cas- 
cade properties in Sec. 14.61 In Sec. \5\ I discuss limitations of the 
model, applicability of physical approximations used in previous 
works, and implication of the results for physics of radiopulsars. 



2 GENERAL NUMERICAL ALGORITHM 

In electromagnetically driven pair cascade in NS magnetosphere 
the following physical processes determine the behavior of the sys- 
tem: 

(i) Charged particles - electrons and positrons - are accelerated 
by the electric field induced by NS rotation. 

(ii) Particles emit high energy gamma-rays. The radiation 
mechanisms relevant for pulsars include curvature radiation, in- 
verse Compton scattering (in both resonant and non-resonant 
regime) of thermal X-ray photons emitted b y the NS, and syn- 
chrotron radia tion of freshly created pairs (e.g. ISturneretal.lll995l : 
IZhang & Har ding 2000). 

(iii) Gamma photons propagate some distance and then create 
electron-positron pairs. In pulsar polar cap the dominatin g process 
is sin gle photon pair creation in the strong magnetic field dSturrockl 
Il97 1). In the outer pulsar magnetosphere the dominant process will 
be photon-photon pair creation either on soft photons emitted by 
the N S (thermal X-rays) or soft photons produced in the cascade 
itself dCheng et al.|[l986h . 

(iv) Creation of electron-positron pairs increase plasma density 
and changes the electric field: if a pair is created in a region with 
strong electric field, electron and positron are accelerated in op- 
posite directions and redistribution of the charge density alters the 
accelerating electric field. 

Probably the best numerical technique for self-consistent 
modeling of plasma kinetics - acceleration of charged particles 
and changes of electromagnetic fields induced by their motion 
(item s 1,4 in the list) - is Particle-In-Cell fe.g. .Birdsall & Langdon 
[1983). There particle distribution is modeled directly by represent- 
ing plasma by an ensemble of macroparticles. PIC is a mature 
numerical technique, many of its properties are well k nown and 
are subject of constant ongoing investigations (e.g. Verb oncoeuJ 
[2005>) . Although on the current stage of the project ID model- 
ing is used, PIC allows straightforward generalization for multi- 
dimension. Particle emission and creation of electron-positrons 
pairs - a radiation transfer problem (items 2 and 3) - in a system 
with strongly and rapidly changing particle energy d istribution are 
best to model utilizing Monte Carlo technique (e.g. ISobol 1 l 19731 : 
iFishmanll 19961) : the computational costs of Monte Carlo are almost 
the same for ID and multidimensional cases. On the other hand, in 
PIC plasma is already represented by discrete particles, what makes 
Monte-Carlo a natural choice. For modeling of pair cascades I de- 
cided to develop a new hybrid PIC/Monte Carlo code. The existing 
hybrid codes used for modeling of gas discharges do not include 
radiation transfer and account only for interaction between charged 



and neutral particles; Monte Carlo technique there is used to ac- 
count for interaction randomness. 

The mean free path of gamma-photons does not depend on 
the plasma density in the polar cap, it is set by the strength of 
the magnetic field and by the curvature of magnetic field lines. 
For the minimum mean free path of photons the estimate of 
iRuderman & Sutherland! (Il975h can be used, which gives A^fp ~ 
10^ cm. In space charge limited flow models photon mean free path 
could be comparable to the NS radius, /Imfp ~ 10^ cm. Characteris- 
tic plasma scales are of the order of the Debye length which depend 
on plasma density. A rough estimate for the Debye length can be 
made assuming plasma density being equal to the GJ number den- 
sity fZoj = T]G}/e: 

,-.4=,(l^!^p.25->/V-'/^cm, (1) 

Ojf \ Me I 

where Bu is the pulsar magnetic field in units of 10^^ G and P is 
the pulsar period in seconds. The photon mean free path is much 
larger than the Debye length, and so it sets the macroscopic scale 
- the length of the computational domain. The Debye length of the 
plasma sets the microscopic scale of computations - the cell size. 
It will be unwise to advance photons in space at the same pace 
as particles - photons propagate large distance to the absorption 
point without interaction while particle motion can change on very 
small spatial scales. Propagation of photons must be calculated sep- 
arately, with larger spatial steps. 

For modeling of electromagnetic cascades in NS magneto- 
spheres I developed a general algorithm which calculates plasma 
motion and photon propagation in diff'erent numerical pace. The 
scheme of the algorithm is presented on Fig.[T] where the sequence 
of operations performed at every time step is shown. 

The plasma dynamics is done with the standard PIC algorithm. 
Using the current density known from the previous step I solve 
Maxwell equations and get electric field at grid points. Then for 
each particle I interpolate the electric field to the particle's position 
and get the electric force on the particle. Solving the equation of 
motion I advance particle momenta and positions. Particle motion 
through the cell boundary is counted as its contribution to the elec- 
tric current. The electric current for each cell boundary is computed 
simultaneously with particle motion and is stored for the next time 
step. 

Photon emission and pair production are calculated as follows. 
I sample how many photons capable of producing electron-positron 
pair each particle emits during the current time step. For each emit- 
ted photon its energy is sampled from the spectral energy distri- 
bution of the corresponding emission process. Then I sample the 
distance the photon will travel until it is absorbed. Calculation of 
the optical depth to pair creation is done with the space steps ad- 
justed according to the current value of the cross-section for photon 
absorption; most of the steps are much larger that the cell size. In 
this way photon propagation is done in the (appropriate) and much 
faster numerical pace than particle advance. Photon's energy, po- 
sition and time of absorption are stored in an array. At every time 
step I iterate over the photon array and pick up photons which are 
absorbed at the current time step. For each of the selected photons 
I inject an electron and a positron at the point of photon absorption 
and delete that photon from the array. Being injected at the same 
point freshly created electron and positron do not contribute to the 
charge and current densities at the time step of injection. 

If there are too many particles of particular kind in the com- 
putational domain, their number can be reduced by deleting some 
randomly selected particles. The total statistical weight of the se- 
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Figure 1. Code structure - sequence of operations performed at every time step. 



lected particles is stored and then statistical weights of all remain- 
ing particles of the same kind are increased in order to compensate 
for the deleted particles. Although this conserves the overall charge 
of the system, the resulting charge distribution will be slightly dif- 
ferent from the one before particle deletion. To proceed with charge 
conserving algorithm one need to solve Poisson equation in order 
to bring the electric field in accordance with the altered charge dis- 
tribution. When the number of photons is reduced, the later step is, 
of course, unnecessary. 



3 ONE-DIMENSIONAL DISCHARGE IN PULSAR 
POLAR CAP 

As previously there were no truly self-consistent studies of electro- 
magnetic cascades - allowing time-dependence and incorporating 
all classes of relevant microscopic processes - 1 decided to address 
first the simplest case in order to develop an intuition about physics 
of pair plasma generation. It was not clear a priori what is the pat- 
tern of plasma flow, and in order to develop an appropriate numer- 
ical technique to model a realistic system with many microscopic 
processes at play a "bare-bone" model must be studied first. 



3.1 Physical model 

The lRuderman & Sutherland! (Il975l) model for pair cascade in the 
polar cap of pulsar is the simplest possible model for a pair cascade. 
Ruderman and Sutherland considered the case when the NS angu- 
lar velocity is anti-parallel to its magnetic momentum - so that the 
Goldreich- Julian charge density is positive - and assumed that the 
work function to extract a positive ion from the surface of a NS is 
much larger than the available electric potential. In this model there 
is no plasma inflow from the surface of the NS and all plasma in the 
cascade zone is produced by pair creation in a series of 'discharges' . 
When enough plasma is produced in the discharge zone, it screens 
the accelerating electric field and, therefore, stops particle acceler- 
ation and pair creation. Plasma flows into the magnetosphere and - 
as there is no source of plasma than (now suppressed) pair creation 
- the plasma density decreases. When there are not enough charged 
particles to screen the accelerating electric field the pair formation 

starts again. 

Although there are hints (e.g. iMedin & Lail 120071) that the 



work function in the NS crust can be small enough - so that par- 
ticles can be extracted from the star - and cascade operates i n the 
so-ca lled space-charge limited flow regime (lArons & Scharle mannI 
|l97i), I think that studying of the Ruderman-Sutherland model is 
worth the eff'ort. As reasons for this I name: i) this model is intrin- 
sically non- stationary and should be a good test of whether non- 
stationary cascade is indeed flexible enough to adjust itself to any 
current density required by the magnetosphere; ii) being the sim- 
plest model it could be used as a testing ground for the numerical 
technique; iii) boundary conditions in this model (no plasma in- 
flow) are similar to that in the problem of electron-po sitron pair 
plasm a generation near the horizon of a black hole jBeskin et al.l 
ll992l : [Hirotani & Okamoto.. 199 8), and, hence, from solution of the 
pulsar problem it would be possible to get some hints to how to 
address the latter problem. 

Ruderman & Sutherlai^ (Il975h estimate the height of the cas- 
cade zone for young pulsars (their eq. (22)) as 

Ks-5x \&p]I^P^'^B\^'^ cm , (2) 

where is the pulsar magnetic field in units of 10^^ G, p6 is the 
radius of magnetic field line curvature normalized to 10^ cm and P 
is the pulsar period in seconds. For young pulsars, with the period 
of the order of ~ 0. 1 sec, h^s is less that the width of the polar cap 

rpc ^ 1.45X lO^p-^/^cm. (3) 

Therefore, ID approximation should work well for such cascades. 
In the Ruderman-Sutherland model the charge density deviate 
strongly from the GJ charge density what creates accelerating elec- 
tric field comparable to the vacuum electric field. The general rel- 
ativistic eff'ects introduce corrections to the electric field of the or- 
der of several per cents of the vacuum electric field (lBeskinlll99d : 
iMuslimov & Tsvganlll990t) , and for this problem they can be ig- 
nored. 

For young pulsars the dominant emission process in terms of 
number of pair-produ ction capable p hotons is the curvature radia- 
tion (e.g. Hibschman & Aronsll2001al) . In this paper I am primarily 
interested in dynamics of the discharge zone, the region with the ac- 
celerating electric field. The size of that zone should be of the order 
of few /zrs, which is of the order of the mean free path of curva- 
ture photons. Synchrotron photons emitted by the injected pairs are 
much less energetic than curvature photons and, therefore, the mean 
free path of synchrotron photons is much larger than h^^. They are 
absorbed at large distances from the NS where plasma density is ex- 
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pected to be very high and electric field is already screened. Hence, 
the pairs produced by the synchrotron photons do not influence the 
discharge dynamics and synchrotron emission can be ignored. 

So, the minimal physical model for the Ruderman-Sutherland 
cascade includes ID electrodynamics, curvature radiation as the 
photon emission process, and pair creation in a strong magnetic 
field as the source of electron-positron pairs. 



3.2 Main equations 

In the superstrong magnetic field of pulsar charged particles are in 
the first Landau level and move strictly along magnetic field lines. 
The radius of curvature of magnetic field lines p is much larger 
than the polar cap radius r^c ; for distances comparable to the width 
of the polar cap particle dynamics can be considered as a motion 
along straight lines. The curvature of the field lines is essential for 
photon emission and pair creation. The radius of curvature of mag- 
netic field lines enters in the expressions for curvature radiation and 
gamma-ray absorption cross-sections as a parameter, which can de- 
pend on particle position. 

I assume that charged particles moves along straight magnetic 
field lines which are perpendicular to the NS surface. A coordi- 
nate axis X is directed along the field lines, its origin is at the NS 
surface and positive direction is toward the magnetosphere. In the 
one-dimensional model charged particles are represented by thin 
sheath with infinite extend in the direction perpendicular to the x- 
axis. I normalize particle momentum to m^c - the normalized parti- 
cle momentum p is its 4- velocity p = j3y, where j3 = v/cis particle 
velocity normalized to the speed of light and 7 = (1 is the 

Lorentz factor. The equation of motion for a particle / is 

dxi 



dt 



(4) 
(5) 



MeC Mi 

qi and m/ are particle charge and mass in units of electron charge e 
and mass correspondingly. W^^ is the term responsible for radia- 
tion reaction. For curvature radiation it is given by 



2e^ / 
3meC rhi 



(6) 



For low energy particles radiation reaction becomes very small and 
for them W^j; in eq. ^ is ignored (see Sec. 13.31 ). 

In one-dimensional model the only changing component of 
electromagnetic fields is the electric field component E parallel to 
the X-axis. The system is essentially electrostatic and the electric 
field can be obtained from the solution of the Poisson equation for 
the electric potential 



— = -47T(T] - T/gj) 



(7) 



SLS E = -dcp/dx. Here t] is the charge density and t/gj is the GJ 
charge density. In order to solve eq. one has to specify either the 
potential diff'erence across the domain or one has to set the electric 
field to some fixed value at one end of the domain, either on the NS 
surface or at the base of the magnetosphere; these boundary condi- 
tions must be specified at every time step. The main free parameter 
in the problem is the average current density which flow trough the 
cascade zone. Charge can accumulate on the NS surface and par- 
ticles can be send back from the magnetosphere. Hence, boundary 
conditions are diff'erent at every time step - the potential drop along 
the computation domain as well as the electric field at the domain 



boundaries change with time. Boundary conditions are related in 
some complicated way to the requirement of providing a certain 
value for the average current density. However, if charge conserva- 
tion is taken into account, the electrostatic boundary value problem 
can be transformed into an initial value problem, which does not 
require boundary conditions. I do so in Appendix lAl where I derive 
the equation for the electric fielcfl 



dE(x, t) 
dt 



-An{j{x, t)-u{t)) 



(8) 



Here j is the current density at point x at time t, and jm is the av- 
erage current which flow through the calculation domain, deter- 
mined by the twist of the magnetic field imposed by the global 
stress balance of the magnetosphere. To solve this equation only an 
initial configuration of the electric field in the domain E(x)1=q is 
necessary; the boundary conditions are incorporated in jm, i.e. the 
electric field at domain boundaries will adjust itself to provide the 
required average current density jm- 

The solution of the eq. ^ gives the correct electric field - the 
one which satisfies the Maxwell equations (in the ID case it is the 
Gauss law) - if one starts from a configuration where E is obtained 
as a solution of the Poisson equation and numerical algorithm con- 
serves electric charge. In my simulations, at the very first time step 
I set some boundary conditions on the electric field (or the poten- 
tial drop in the domain) and some initial particle distribution; then I 
compute the charge density rj and solve the Poisson equation ^ for 
that boundary conditions to get the initial distribution of the electric 
field in the domain. At all subsequent time steps for each point in 
the numerical grid I compute the electric field from its value at the 
previous time step using equation (|8l); the current density j due to 
particle motion is calculated using a change conserving algorithm. 

Physically, in order for electric field to be zero in the polar cap 
the charge density must be equal to the GJ charge density and the 
current density equal to the jm required by the magnetosphere. The 
GJ charge density enters in the Poisson equation which is solved 
at the first time step; because of charge conservation the system 
"keeps memory" of the GJ charge density at all subsequent time 
steps. The current density enters in the equation for the electric 
field. So, the system tries to adjust to both these requirements. 

Particles moving along curved magnetic field lines emit pho- 
tons via curvature radiation mechanism. Spectral energy distribu- 
tion of curvature photons emitted by a particle with th e Lorentz 
factor 7 is given by the standard formula (|jacksonlll975h 



^A^ph^^^_ 1 cufC 
dtde V3;r Xc 7 



1 r 



(9) 



where 6 is the photon energy normalized to m^c^, af is the fine 
structure constant, Zq = h/nieC = 3.86 X 10~^^ cm is the reduced 
Compton wavelength, ^5/3 is the modified Bessel function of order 
5/3; e^^^ = (3/2)/lcP~V^ - 57.92 Pg the peak energy of cur- 
vature photons, 76 = 7/10^. The integral in eq. ^ has asymptotic 
forms 



15^-2/3 - 1.81, ify<^ 1 
25^-V^^^ ifj»l 



(10) 



The total number of curvature photons with energies greater that 
some Sa emitted by the particle during time dt is 



^ iLevinson et al.l bOOSh used the same approach for calculation of the elec- 
tric field. They did not elaborate on the physical meaning of j^, so I decided 
to present a detailed derivation of eq. (8} 
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1 (XfC 1 / 6 

dN^^{e> ea)^dt — ^ -F\- 

where 



(11) 



(12) 



dxKs/3(x). 

For small values of its argument F(f ) has the following asymptotic 
form 



F(0 ^ 1 + 0.346^ - f^/^(1.232 + 0.033^^), ^ « 1 . 



(13) 



Only very high energy photons capable to produce an electron- 
positron pair in the calculation domain are of relevance for the 
considered problem and only those are tracked in the code (see 
Sec.ES. 

I assume that photons are emitted tangentially to the magnetic 
field lines and then move along straight lines. The angle if/ between 
the photon momentum and the magnetic field increases as the pho- 
ton propagates further from the emission point. In a simple model 
where magnetic field lines have constant curvature the angle be- 
tween the photon momentum and the magnetic field is given by 



1//(X) = (X-Xe)/P, 



(14) 



where Xe is the coordinate of the emission point. In the dipolar mag- 
netic field the expression for i//(x) is slightly more complicated 



3 X - ^ 

4 X 



1 + ■ 



(15) 



here 6e is the colatitude of the emission point (a free parameter in 
the ID model), R^s is the NS radius. Both models were used in 
the singulatio ns. The cross-section of photon absorption is given by 
(lErbei]ll966h 



af B , ^ 

o-By = 0.23 — — smj/r exp - — 



(16) 



where = eBsmij/IBq ^ 2.27 X lO-^e^is^A, and B^ = m]c'len ^ 
4.41 X 10^^ G is the critical magnetic field strength. The cross- 
section grows exponentially as photon propagates further from the 
emission point. 

When the photon is absorbed I assume that its energy is 
equally divided between newly created electron and positron. The 
perpendicular component of particle's momentum will be rapidly 
radiated as synchrotron photons, which, as described before, are 
neglected. Injected particle ends up having only the longitudinal 
component of the momentum 



1/2 



(17) 



where i//a is the angle between the photon momentum and the mag- 
netic field at the absorption point. 



3.3 Numerical implementation 

In this section I describe normalization of physical quantities, intro- 
duce several numerical parameters controlling the algorithms, and 
give their typical values in my simulations. I also give an overview 
of main numerical algorithms used in the code; a detailed descrip- 
tion of the numerical code will be given elsewhere. 

Distances are normalized to the radius of the pulsar polar cap, 
•^0 = r^c^ given by eq. The electric potential is normalized to 
the vacuum potential drop between the rotation axis and the edge 
of the polar cap in the aligned rotator 



Oo = - ^ ^ 6.6 X 10^^ BuP~^ V , 
c 2 



(18) 



Q is the pulsar angular velocity. The electric field is normalized to 



(19) 



^0 = ^ ^ 4.6 X \O^BnP~^^^ V/cm, 

Xo 

and the charge density to the absolute value of the Goldreich- Julian 
charge density 



^0 



QB 

Ittc 



(20) 



Each numerical particle is a macroparticle representing a 
(large) number of real particles A^o. either electrons or positrons. 
Each numerical particle has a statistical weight w/ = w/A^o- When 
a macroparticle emits a photon, the latter gets the particle's statis- 
tical weight (also see below the description of photon sampling); 
when the photon is absorbed the injected electron and positron get 
the photon's statistical weight. An important numerical parameter 
is N^f^ - the number of macroparticles with the normalized statis- 
tical weight Wi = 1 in a cell which create Goldreich- Julian charge 
density 



^0 



(21) 



The parameter controlling the number of numerical particles in the 
simulation is N^f^; Nq is computed at the start of the simulation 
from eq. (|2TI ). The diff'erence in the number density between parti- 
cles of opposite signs of the order of A^^^^^ results in a large electric 
field; this number should be not very small, otherwise numerical 
noise will strongly contaminate results. In my simulations values 
of N^f^ > 5 provide acceptable level of numerical noise which al- 
lows to recognize plasma oscillation excited in the pair plasma. 

The calculation domain is divided in equal numerical cells; 
a typical value of in my simulations is several thousands. I 
use ID vers ion of the charge conservative algorithm proposed by 
IVillaseno7 & Buneman ( 1992) for scattering of charge and current 
densities to the grid points and for interpolation of the electric field 
to particle' locations. Integration in time of eq. (|4]), (|5]), (|8l) is done 
with a leap-frog scheme with a uniform time step A^. The radia- 
tion reaction in eq. ^ is taken into account only for particles with 
momentum larger than a certain value p^^^. 

For each particle if its momentum is larger than a certain value 
I calculate the mean number of photons A^ph with energies 
larger than a certain €^ the particle emits during time At accord- 
ing to eq. fTTI ). If A^ph is small, less than a certain A^^^^, I sample the 
number of actually emitted photons from the uniform distribution 
with the mean value equal to A^ph. Then for each photon I sample 
its energy from the distribution given by eq. iT2\ using either cut- 
point or inverse transform methods (Fishman 1996). The values of 
F(0 in eq. ([12]) are tabulated for 0.01 < ^ < 10 for use in cut- 
point method, for smaller f invers transform method is used with 
the asymptotic formula F(0 ^ 1 - 1.232^^/3 If A^ph > N^^"", the 
particle emits a fixed number of numerical photons - the spectrum 
is divided in A^^^" bins, the number of emitted numerical photons is 
equal to N^^; each photon gets a statistical weight equal to the prod- 
uct of the statistical weight of the emitting particle and the number 
of photons emitted in the corresponding spectral energy bin given 
by eq. (|9]). 

To calculate the position where the photon is absorbed I sam- 
ple the optical depth the phonon should achieve before being ab- 
sorbed; then I integrate the cross-section ([T6] | along phonon's tra- 
jectory until the required optical depth is reached. The cross-section 
of photon absorption in the polar cap grows exponentially with the 
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distance from the emission point, and most of the trajectory do not 
make significant contribution to the optical depth. At first optical 
depth along photon's trajectory is calculated using rectangle meth- 
ods with large spatial steps (~ 1/20- 1/40 of the domain size) until 
the optical depth on the next step would exceed the required value. 
This integration method overestimates the optical depth, the trajec- 
tory always continues beyond this intermediate stop point. I redo 
the cross-section integration between the emission and the stop 
points using 15 points Gauss-Kronrod integrator what provides a 
very accurate value for the optical depth at the stop point. Then the 
optical depth integration is proceeded with a smaller spatial step 
comparable to the cell size. The number of cross-section evalua- 
tion in this algorithm is, on average, by a factor of few tens smaller 
than a cross-section integration with the step equal to the cell size 
would require. 

Values of the numerical parameters ^ p^^, eZ"", N^^"", A^cr 
are fixed at the start of the simulation. These values are chosen to 
sample all pair creation capable photons and correctly account for 
the radiation reaction on one side, and to minimize the computation 
time of the other side; their particular values depend on physical 
conditions: the pulsar period, the magnetic field strength, and the 
radius of curvature of magnetic field lines. The typical values I used 
in the simulations: p^"" ~ p^^ ~ 10^ - 10^ eZ"" -2-40, 7V^^^^ ~ 
A^cif - 50. 

The numerical code was developed from scratch and written in 
C++ programming language. Its modular object-oriented structure 
is designed to facilitate further extension to multidimension and in- 
corporation of additional physical processes. I tested the PIC part 
of the code performing simulations of the following test problems: 
oscillation of two test particles, two stream instability in both rela- 
tivistic and non-relativistic regime, non-relativistic and relativistic 
Child's laws, dependence of plasm a frequency on numerical res- 
olution teirdsall & Langdonlll985h . I also tested that the code is 
indeed change conservative up to machine precision. The Monte- 
Carlo part of the code was tested as follows. I verified that the en- 
ergy distribution of emitted photons agrees with the spectrum of 
curvature radiation. For several fixed emission points, diff'erent val- 
ues of photon energy, magnetic field strength, and radius of curva- 
ture of magnetic field lines I compared the distribution of photon 
absorption points produced by the Monte-Carlo code with the cor- 
responding theoretical distributions. I also checked that for a given 
time interval the total energy of emitted photons is equal to the par- 
ticle radiation reaction losses. 



4 RESULTS OF NUMERICAL MODELING 

The numerical simulations have shown that in the Ruderman- 
Sutherland model pair creation is quasi-periodic and self- sustained. 
I performed simulations for diff'erent initial particle distributions, 
initial electric fields, strengths of the magnetic field, radii of cur- 
vature of magnetic field lines, and pulsar periods. Independent on 
the initial configuration for non-zero pair-creation process be- 
gins some time after the start of simulations. How and how much 
plasma is formed in this initial burst depends on the specific setup. 
Plasma generation stops after enough plasma is produced to screen 
the electric field. After plasma generated in this burst of pair forma- 
tion leaves the domain (in a couple of domain flyby times), behav- 
ior of the cascade for given magnetic field, pulsar period, and the 
mean current jm is the same, independent on the initial configura- 
tion. All subsequent bursts of pair formation do not depend on the 
initial setup - the system seems to forget the initial conditions. Af- 



ter that initial burst of pair creation the cascade zone always settles 
down to a quasi-periodic behavior. For a given cascade behavior 
is qualitatively similar for all other physical parameters admitting 
pair creation. 

I describe here main properties of the Ruderman-Sutherland 
cascade using as an example a pulsar with period P = 0.2 s, mag- 
netic field B = 10^^ G, and the radius of curvature of magnetic field 
lines p = 10^ cm. The radius of curvature of magnetic field lines 
comparable to the NS radius implies that there is a non-dipolar 
component of the magnetic field in the polar cap region. I per- 
formed simulations for pure dipole magnetic field with p ~ 10^ cm 
too. Qualitatively results do not depend on the radius of curvature, 
but for smaller p calculations with the same numerical resolution 
can be done faster because the size of the gap with accelerating 
electric field is smaller. On the other hand, adoption of this value for 
p simplifies comparison with the original Ruderman & Sutherland! 
(1975) model, where the same value for p was used. 

The polar cap radius for such pulsar is r^c = 3.24 X 10"^ cm 
(eq. JS])). The heights of the gap should be (see eq. (|2])) 

h^s,i ^ 2.5 X 10^ cm ^ O.OUrp, , (22) 

and the potential drop in the gap is 

AVrs,! = 2;r^Gj/^Rs,i ^ 1-98 x lO^^y = 0.0120o , (23) 
so the maximum Lorentz factor of electrons and positrons is 

MVrc 1 

^max ^ ^ 3 ^ ^ (24) 

rUeC^ 

The angular velocity of NS rotation is anti-parallel to the mag- 
netic moment of the star and the Goldreich- Julian change density is 
positive. The length of the computation domain for the simulations 
described in this section is L = 0.3 r^c - 9.72 X 10^ cm. Numerical 
grid has = 5000 points, so that the cell size Ax ^ 1.94 cm. 
The number of numerical particles in cell providing the GJ charge 
density N^f^ = 10. Other numerical parameters are p^^^ = p^^ = 
5 X 10^ eZ"" = 20, N^^"" = 50, N^^ = 80. 

I describe properties of cascade with physical parameters 
given above for 3 diff'erent current densities: = Jgj, jm = 0.5 Jgj, 
and Jm = 1.57gj. First I describe main properties of cascade with 
jm = 7gj- Pair formation dynamics for diff'erent current densities is 
qualitatively similar. Later in this section I will highlight the diff'er- 
ences in cascade properties for = 0.5 joj and jm = 1.5 7gj. 

4.1 Pattern of plasma flow 

In this subsection I describe the pattern of plasma flow for a typical 
cycle of pair formation in cascade with j^^ = 7gj. Cascade develop- 
ment is illustrated by a series of snapshots at several time moments 
during a cycle of pair formation taken from a long simulation where 
several such cycles were observecj^ In Fig.|2]l plot the change den- 
sity at equally spaced time interval during the discharge cycle. In 
the upper panel of that figure I present an overview of the entire 
cycle, in the lower panel I plot snapshots of the change density dis- 
tribution at smaller time intervals for the most interesting part of the 
discharge - formation of a new plasma blob. In Figs.[3l|4]more de- 
tailed information about physical conditions in the discharge zone 

In a previous short publication (lTimokhinll2009l) I presented plots simi- 
lar to Figs. 12] S of this paper for a different cycle of the same simulation. 
Comparing these plots one can see that different bursts of pair formation are 
indeed very similar. 
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Figure 2. Snapshots of charge density distribution in the calculation domain for cascade with jm = Jgj. Charge density 77 as a function of distance x from 
the NS is plotted at equally separated moments of time; 77 is normalized to the Goldreich- Julian change density t/gj. The time t shown in small square boxes 
is normalized to the flyby time of the computation domain and is counted from the start of the simulation. The presented cycle is taken from the middle of a 
long simulation, top: The whole cycle of cascade development, bottom: Snapshots for time interval marked by the gray area in the top panel; these snapshots 
illustrate formation and propagation of plasma blob in more detail. 



is shown: the number densities of electrons and positrons 7/+, the 
accelerating electric field phase portraits (p-x diagrams) of elec- 
trons, positrons and pair producing photons. In the phase portraits 
particles with positive values of 4-momentum p are those which 
move from the NS, particles with negative p move toward the NS. 
The time t in this figures is normalized to the flyby time of the com- 



putational domain - the time a relativistic particle needs to cross the 
domain L/c. The time is counted from the start of a particular sim- 
ulation, so its absolute value has no physical meaning - only time 
intervals between the shots have physical meaning. 

Each pair creation cycle could be conveniently divided into 
three phases: i) vacuum gap formation (timeshots for t = 6.033 - 
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Figure 3. Ignition of pair formation in cascade with jm = ;gj. Several physical quantities are shown as functions of the distance x from the NS. Plots in each 
column (for the same time are aligned - they share the same values of x. Snapshots are take at time moments of the first three marked snapshots in the bottom 
panel of Fig.|2] The following quantities are plotted: row: 77+ - charge density of electrons (negative values, blue line) and positrons (positive values, red 
line); 77+ is normalized to the Goldreich- Julian charge density t/gj. 2"^ row: the total charge density 77 normalized to the Goldreich- Julian charge density t/gj. 
yd YQy^. accelerating electric field E normalized to the vacuum electric field Eq. 4^^ row: phase space portrait of positrons (horizontal axis - positron position 
X, vertical axis - positron momentum p+ normalized to m^c). 5^^ row: phase space portrait of electrons (horizontal axis - electron position x, vertical axis - 
electron momentum p- normalized to m^c). 6^^ row: phase space portrait of pair-producing photons (horizontal axis - photon position x, vertical axis - photon 
momentum py normalized to m^c). 



6.633 in Fig. [3, ii) formation and propagation of a plasma blob 
(t = 6.633 - 7.833) iii) relaxation (t = 7.833 - 8.833). Each burst 
of pair formation generates dense electron-positron plasma which 
screens the electric field. Particles must leave the domain in order 
to provide the required current density. When plasma leaves the po- 
lar cap a gap with almost no particles inside is formed; the vacuum 
electric field in the gap is no longer screened (phase (i)). The few 



particles in the gap are accelerated and emit high energy gamma- 
photons, and the process of plasma creation starts again. Electron- 
positron plasma is produced non-uniformly, it forms a blot0 of rel- 
ativistic plasma where large amplitude plasma oscillations are ex- 



^ Actually I am computing plasma sheets in ID, but in 2D and 3D these 
would be plasma "blobs", so I use the latter term throughout the paper. 
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Figure 4. Screening of the electric field in cascade with jm = Jgj. Snapshots are take at time moments of the last four marked snapshots in the bottom panel 
of Fig. 12] The same quantities are plotted as in Fig. [3] 



cited (phase (ii)); the blob is visible in Fig. [2 as a packet of large 
amplitude charge density oscillations. The blob moves into the 
magnetosphere leaving a tail of moderately relativistic plasma be- 
hind. When the blob leaves the computational domain, the remain- 
ing plasma still screens the vacuum electric field till the plasma 
density drops below na and pair formation starts again (phase (iii)). 
Below I describe these processes in more details. 

A typical cycle starts with formation of a vacuum gap above 
the NS surface (timeshots at t = 6.033 - 6.650 in Figs. E] [3]). The 
gap forms because plasma leaves the domain in order to provide 
the required current density jm- The GJ charge density is posi- 
tive, and when the electric field is completely screened there are 
more positrons than electrons. The current density is positive, and 



positrons, on average, must flow toward the magnetosphere (to the 
right). As there is no plasma inflow through the domain bound- 
aries, the gap forms when there are not enough charged particles 
to provide the required current and change densities in the whole 
domain; above the gap the plasma still sustain the required current 
and change densities. Positrons flow into the magnetosphere, so the 
gap forms at the NS surface. 

Pair creation starts close to the NS and is ignited by gamma- 
rays emitted by electrons flowing toward the NS. These primary 
electrons have been created in the previous burst of pair formation, 
they leak from the tail of the plasma blob formed in the previous 
cycle and enters the gap from above. These electrons are visible 
as a thin line of particles with negative momenta /?_ in the elec- 
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tron phase space portraits (the 5^^^ row in Fig. [Sj. Pair production 
capable gamma-photons emitted by these electrons have negative 
momenta and are visible in the photon p - x diagrams as scattered 
dots with negative Py at t = 6.5,6.6. These photons are eventually 
absorbed and create electron-positron pairs, which are visible as 
scattered dots at the left in the electron and positron p-x diagrams 
at ^ = 6.6 - 6.7. These newly created electrons and positrons are 
accelerated by the strong electric field of the gap (see timeshot at 
t = 6.6); when they have been accelerated up to Lorentz factors of 
the order of several 10^ they start to emit pair production capable 
photons. Att = 6.700 there are already gamma-photons with pos- 
itive momenta, they have been emitted by the secondary positrons 
accelerated in the strong vacuum electric field. At these early stages 
of pair discharge the density of the newly created plasma is still 
very low (see plots for 77+, the row in Fig. [3]), and the electric 
field is not influenced by the injected plasma (the 3^ row in Fig. [3]). 
All electrons and positrons are accelerated up to high Lorentz fac- 
tors and start emitting pair production capable gamma-photons very 
soon. 

The secondary electrons and positrons are accelerated in op- 
posite directions, and plasma start being polarized (see distribution 
of the charge density at ^ = 6.8, 6.85, 6.9 in Fig.[2]). The polarization 
of plasma creates an electric field opposite to the vacuum field, and 
the eff'ective accelerating electric field decreases. When the particle 
number density become comparable to the Goldreich- Julian density 
the accelerating electric field starts being screened by plasma. 
How the electric field is screened is shown as a series of snapshots 
in Fig. m The screening naturally starts in the place where plasma 
density is maximal. When more and more plasma is injected the 
region of screened electric field broadens till it eventually extends 
up to the NS surface. 

The particles which produce the most of the pair creating pho- 
tons are the secondary positrons which have been accelerated at the 
time when plasma density was small and electric field was strong. 
Secondary electrons have been accelerated up to very high ener- 
gies too, but they have been moving toward the NS - they slammed 
into the NS surface and do not contribute to pair creation at later 
times. The high energy positrons move into the magnetosphere 
emitting gamma-rays which turn into electron-positron pairs. These 
positrons practically co-move with the gamma-rays. Freshly cre- 
ated pairs are relativistic and have momenta directed from the NS; 
most of them will remain relativistic and will move into the magne- 
tosphere. So, the pair plasma forms a blob with constantly increas- 
ing particle density (see the 1*^ row in Fig.|4l). 

Photons cannot go ahead of relativistic particles at the front 
edge of the blob, and pairs are injected only inside the blob. As 
the plasma blob moves into the magnetosphere so does the vacuum 
gap limited from below by the front edge of the blob. Ahead of the 
blob there is practically vacuum with strong electric field, as there 
are only electrons leaking from the tail of the previous blob and 
their number density is very low. The electric field inside the blob 
is screened by the plasma. At the front of the blob a sheath of posi- 
tive charge is formed which screens the vacuum electric field. This 
sheath is visible as a large spike in the charge density distribution 
at timeshots t = 7.05 - 7.433. Inside the sheath the electric field 
goes from the vacuum value to its very low value in the blob. Pairs 
injected in this sheath by conversion of gamma-photons emitted by 
the particles which are already in the sheath are accelerated by that 
electric field and start emitting gamma-rays too. As more and more 
pairs are ejected there the width of the sheath decreases. However, 
the number of particles in the sheath is small (see plot for 77+ in 
Fig- 111 at t = 7.25). Pairs injected in other parts of the blob are not 



accelerated because the electric field there is screened. Hence, sub- 
sequent pair formation is driven mostly by particles accelerated at 
early stages of the discharge, when the electric field was strong. As 
these particles are emitting gamma-rays and they are not acceler- 
ated anymore, their kinetic energy decreases. This can be seen in 
the positron p-x diagrams in Fig.|4j the spike of high energy par- 
ticles at the blob's front is due to particle acceleration in the charge 
sheath. 

I observe thermalization of freshly created electrons and 
positrons in the simulations. Low energy particles are present start- 
ing from very early stages of blob formation. Some of the low 
energy particles are reversed back. While the bulk of the plasma 
moves with relativistic speed into the magnetosphere some parti- 
cles are left behind forming a "tail" of the blob which has much 
lower particle number density than the blob itself. Although the 
fraction of particles left behind is small, their number is enough to 
screen the vacuum electric field for some time, preventing imme- 
diate formation of a new vacuum gap after the blob detaches from 
the NS. 

While the general structure of the flow is evident from the per- 
formed simulations, some questions remain unanswered. The most 
important among them is about the time between discharges. It de- 
pends on the rate of plasma leakage from the blob - the more par- 
ticles leak out, the later the next gap forms. Due to continuous pair 
injection plasma density in the blow increases enormously and at 
some time the numerical scheme stops resolving the Debye length 
of the plasma, and results start depending on the numerical resolu- 
tion. Because of this the blob cannot be followed for time interval 
long enough (and distances large enough) to get the repetition rate 
of the cascade. In the presented simulations the size of the simula- 
tion zone is set such that the blob leaves calculation domain before 
the numerical scheme fails to model it correctly. When the blob 
leaves the calculation domain particles are still leaking from it into 
the tail. When the blob is no longer in the computational domain, 
particle supply to the tail is stopped and the time interval during 
which plasma density in the domain drops to the GJ density - and a 
new gap begins to form - is substantially smaller in my simulations 
that it would be in reality. Duration of the relaxation phase in the 
simulations (timeshots t = 7.833 - 8.833) is strongly influenced by 
the numerical setup. 

However, I believe that the qualitative behavior of the plasma 
during this phase is represented correctly, as there are no physi- 
cal processes in the blob tail except the particle supply from the 
blob that the code fails to model when the blob is no longer in 
the computational domain. There is no particle acceleration in the 
tail and pair creation is only due to electrons which leak from the 
previous blob and are accelerated toward the NS in the traveling 
vacuum gap. The number of these electrons is negligible compared 
to the number of pair-producing positrons in the blob. The fraction 
of particles leaking from the blob when it is still inside the domain 
is also small, of the order of few per cents. This is enough to screen 
the vacuum electric field, but the energy carried by those particles 
is negligible compared to the energy carried by particles in the blob 
(see Sec. 14.41 ). Hence, the only cascade characteristic substantially 
influenced by the rate at which particles leak from the blob is the 
repetition rate of pair creation bursts. Already from this simulations 
it is clear that the time between discharges is longer than the vac- 
uum gap crossing time h^s/c. This introduces a new time scale into 
the Ruderman-Sutherland model. 

Qualitatively the plasma flow after the discharge should be as 
follows. The tail consists of particles leaked from the blob; those 
are mildly relativistic particles, some of them are trapped in plasma 
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oscillations, and, on average, the tail moves with a subrelativistic 
velocity. The vacuum gap is limited from above by the tail of the 
previous blob: the plasma must move in order to support the current 
density when the plasma density drops to values comparable to 
the GJ density, the plasma cannot support the required current den- 
sity and the GJ charge density at the same time, a gap appears, and 
so the tail ends in a vacuum gap. Among the trapped particles in the 
tail there are both electrons and positrons which move toward the 
NS. The electrons at the tail's end enter the gap and get accelerated 
toward the NS - they will be the primary particles in the next burst 
of pair formation. The positrons entering the gap are reversed and 
are sent into the magnetosphere. The upper boundary of the vac- 
uum gap - the tail's end - moves with a subrelativistic velocity, the 
front of the blob is formed by ultra-relativistic positrons and moves 
relativistically; therefore, the gap shrinks when the blob moves into 
the magnetosphere. The velocity of the blob's tail Vtaii is several per 
cents less than the speed of light, Vtaii ~ 0.95 -0.99c. Eventually the 
front of the new blob catches the tail of the previous blob; the dis- 
appearance of the gap is directly visible in the case of - 0.5 Jgj 
(see Sec. 14.51 ). Therefore, the magnetosphere in the open field line 
zone will be filled with plasma everywhere and starting at some 
distance from the polar cap there will be no gaps in plasma spacial 
distribution. 

4.2 Superluminal plasma wave 

When plasma starts being injected into a region with strong electric 
field it is polarized and starts screening the electric field. During the 
process of screening of the vacuum electric field large- amplitude 
oscillations are excited in the injected pair plasma; these oscilla- 
tions are visible in timeshots starting dXt - 6.85 and until the blob 
leaves the domain. Screening of the electric field starts in the mid- 
dle of the blob and spreads to its edges. This spreading occurs in 
the form of an electrostatic wave. The propagation of the wave can 
be seen in Fig. [5] where I plot snapshots for the electric field, the 
charge density, and the particle number density for the same spatial 
domain where the blob is being formed for 6 moments of time; I 
plot also three vertical lines which mark fiducial positions moving 
with the speed of light toward the magnetosphere. One can clearly 
see that the phase velocity of the wave is greater that the speed 
of light; also note, that the wave propagating toward the NS is su- 
perluminal too. In the process of electric field screening less and 
less charge separation would be necessary to kill the electric field. 
Apparently this is the reason why the wavelength of plasma oscilla- 
tions decreases. At the time of wave formation the particle number 
density is already very high, there are more than ~ 100 numerical 
particles per cell. The Debye length of the plasma - calculated as 




where n{y) is the number density of particles with the Lorentz 
factor 7 - is resolved; at the time the snapshots shown in Fig. [5] 
are made /Id is several tens times larger than the cell size. Hence, 
these oscillations are not numerical artifacts. The electric field in 
the wave is too weak to accelerate particles up to energies when 
they can emit pair production capable photons. Particle injection 
rate is set by very high energy particles accelerated at earlier time 
and there is no back reaction of the wave on the particle production 
rate. 

Decreasing of the wavelength eventually ends when the wave- 
length become equal to the cell size. From that moments the code 
cannot correctly follow propagation of the wave. In timeshots for 



^ > 7.25 the wavelength is not resolved anymore. Oscillations per- 
sists, but because the wavelength is equal to the cell size, param- 
eters of the wave depend on the numerical resolution and so tell 
us little about how the wave would propagate at that time in a real 
cascade. However, based on the information obtained at the early 
stages of wave evolution (for t < 7.25) - when results of numerical 
modeling should be reliable - I would like to make the following 
remarks. Being superluminal this wave is not damped via Landau 
damping. I may speculate that this wave could stay superluminal, 
with its phase speed approaching the speed of light from above, 
for a time long enough that it can travel into the magnetosphere 
practically undamped. Although in my ID simulations the wave is 
electrostatic, in reality it would be electromagnetic. It could escape 
the magnetosphere and be observed as coherent pulsar emission, 
or/and it can excite another electromagnetic wave(s) which escape 
the magnetosphere. 

4.3 Particle momentum redistribution and current 
adjustment 

Starting from very early stages of blob formation the injected pairs 
start being thermalized. I use the term "thermalization" here rather 
loosely, meaning that in the particle momentum distribution there is 
a strong broad component which peaks at some momentum value; 
it extends up to very low energies decreasing like a power law, and 
decreases strongly after the peak. Although I did not perform a for- 
mal fitting procedure for particle momentum distributiorQ, the low 
energy tail of the "thermal" component follows the ID Maxwell- 
Juttner distribution dn/dp ~ const for small p quite good. In Fig.[6l 
I plot particle momentum distribution p (dn/dp) for three diff'erent 
moments of time. In the upper panel I plot the momentum distribu- 
tion of particles moving toward the magnetosphere, p is positive; in 
the lower panel - the momentum distribution of particles moving 
toward the NS, /? is negative. These distributions are for particles 
in the blob - they are averages over x e [0,0.135] for t = 6.95, 
X G [0.05,0.18] for ^ = 7.1 and x e [0.1,0.225] for t = 7.25 (cf. 
plots for the particle number density in Fig.|4l). There are low en- 
ergy particles and there are particles with both direction of motion 
in the blob. Essentially the pair plasma become four-component: 
(i) positrons moving into the magnetosphere, (ii) positrons moving 
toward the NS, (ii) electrons moving into the magnetosphere, (iv) 
electrons moving toward the NS. Such four-component plasma eas- 
ily adjusts locally to both requirements of providing the GJ charge 
density and the current density j^. 

At initial phases of blob formation at least one of the pro- 
cesses leading to plasma thermalization is trapping of particles in 
the electric field of the plasma wave. In Fig.[7]trajectories for three 
of such trapped pairs are shown. In that figure phase trajectories of 
pair-producing photons are plotted by dotted lines, marked as 71,2 3- 
When photon is absorbed an electron and a positron are injecteqj. 
Electron trajectories are shown by dashed lines, positron trajecto- 
ries by the solid lines; trajectories' final points are marked as ^^23 

^ To get an accurate momentum distribution it is necessary to collect 
enough numerical particles. This leads to averaging over a macroscopic vol- 
ume with different physical conditions, and results of such fitting would be 
ambiguous anyway 

^ Note that the initial kinetic energy of injected pairs is much less than 
the energy of the pair producing photons; this is easy to see although from 
eq. fTTl . The rest of the energy goes into synchrotron radiation; that energy 
is carried by many photons with energies much less than that of the primary 
photon. 
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Figure 5. Formation and propagation of superluminal electrostatic wave in the forming plasma blob for cascade with jm = Jgj • There are six snapshots for 
the electric field E, the total change density 77 and the charge density of electrons (negative values, blue line) and positrons (positive values, red line) 77+. All 
quantities are plotted as functions of distance x for the part of the calculation domain where the blob is forming. Snapshots are taken at equally separated time 
intervals. Plots in each column are aligned and share the same values of x. The same normalizations for physical quantities are used as in Fig. [3] The three 
thin red vertical lines in each plot mark fiducial points moving with the speed of light toward the magnetosphere. The wave is superluminal, its maxima move 
faster that these lines. 



and ^[23 correspondingly. Particle trajectories end at t = 7.1; at 
that and earlier time both the Debye length and the wavelength of 
the wave are well resolved, they are tens times larger than the cell 
size, and there are many particles per cell; particle trajectories are 
well resolved too. Hence, the thermalization is not a numerical arti- 
fact. Thermalization of freshly injected pairs proceeds also at later 
stages of cascade development, however, as the code does not re- 
solve the plasma wave anymore, it is not possible to disentangle the 
influence of the wave on the thermalization process at this time. 

From the current simulations it is not clear what is the fate 
of the plasma wave in the blob. If it exists for a long time, parti- 
cle thermalization in these oscillations could continue. On the other 
hand, deviation of the current density from results in appearance 
of an electric field even if the charge density is equal to the local 
GJ charge density. The presence of that electric field in a plasma 
with broad momentum distribution might result in some instability 
which could facilitate pair thermalization - at least I see plasma 



thermalization during the relaxation phase when there is no large 
amplitude plasma wave and plasma leaves the domain adjusting to 
the required current density by redistribution of particle momenta. 
This topic, however, needs additional investigation and will be ad- 
dressed in future publications. 

To reverse the direction of motion of low energy particles a 
weak electric field would be sufficient. For charge density to be 
equal to the local GJ charge density the number of positrons should 
be greater that the number of electrons by Hq}. In order to provide 
current density less than Jgj some positrons should move toward 
the NS. If there is a population of low energy trapped positrons, 
a weak fluctuating field can ensure that some trapped positrons on 
average will be moving toward the NS. A similar process can pro- 
vide a current density larger than Jgj and still keep the change den- 
sity equal to the GJ charge density. The adjustment of the current 
density can proceed in the cascade zone locally, without inflow of 
charged particles from the magnetosphere (the latter mechanism of 
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Figure 6. Particle momentum distribution for plasma in the blob at three moments of time for cascade with jm = Jgj • Positron distributions are shown by solid 
lines, electron distributions - by dashed lines, distribution of pair producing photons - by dotted lines. Plots in the top row show distributions for particles 
moving toward the magnetosphere {p > 0), plots in the bottom row - distributions for particles moving toward the NS {p < 0). Each column corresponds to 
the same moment of time shown above the plots. Plots in each columns are aligned and share the same values of \p\. The following blob sizes were assumed: 
X e [0, 0.135] for t = 6.95, x e [0.05, 0.18] for f = 7.1 and x e [0.1, 0.225] for t = 7.25. 



curren t adjustment was discussed in iLvubarskijl ( ll992h : lTimoldiiiil 

In Fig. [8] I plot electric currents through the lower domain 
boundary toward the NS surface (dashed line, this current should 
be negative) and through the upper end of the calculation domain 
(solid line, this current should be positive) as functions of time. The 
required current density jm is achieved on average, at each moment 
of time the current density deviates from jm. Fluctuations are the 
strongest when particles from the burst of pair formation hit the NS 
surface, and when the blob reaches the outer boundary. In all cases 
the relative deviation of the mean over the cycle current density 
from jm is less than ~ 10~^. 

4.4 Cascade energetics 

The height of the gap is ~ 2 times larger than the estimate given 
by eq. (|22] ). The reason for this is that the upper boundary of the 
vacuum gap - the end of the blob tail - is moving into the magne- 
tosphere while the electrons which ignite the cascade are moving 
toward the NS. When these electrons arrive at the point where they 
emit first pair production capable photons, the gap's upper bound- 
ary has moved some distance into the magnetosphere and the gap 
size is larger. 

The electric field in the gap linearly increases toward the NS 
and the first secondary particles are injected into the region with 
very strong electric field. A substantial amount of particles needs 
to be generated before the vacuum electric field is screened. In 
the meanwhile the vacuum gap is still growing and freshly in- 
jected particles are accelerated in a very strong electric field. In 
the considered case a noticeable number of particles reaches the 
radiation-reaction limited Lorentz factor ~ 1.4 X 10^ (see snapshots 
at t = 6.8, 6.95 in Fig.g]). This energy is ~ 4 times higher than y^^'l 
given by eq. (l24l ). When the electric field is screened these particles 
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Figure 7. Trajectories of three pairs and their parent photons in the phase 
space (distance x is along the horizontal axis, momentum p - along the ver- 
tical axis) for cascade with = ;gj- Trajectories of pair creating photons 
are shown by dotted lines and marked as 71,2,3- Trajectories of positrons 
are shown by solid lines and marked as trajectories of electrons are 
shown by dashed lines and marked as 3- Marks are located at the ends 
of corresponding particle trajectories. All trajectories end at f = 7.1. 



start losing their energy quickly and then for the most of the pair- 
producing positrons the Lorentz factor do not exceed ~ 8 x 10^ 
(see snapshots at ^ = 7.1, 7.25 in Figs.|4]). yi is the Lorentz factor of 
a relativistic electron/positron which loose substantial amount of its 
kinetic energy due to curvature radiation while moving a distance 
comparable to the length of the computation domain: yi/^t ~ W^', 
6t - L/c and W^r is the radiation-reaction term in the equation of 
motion (|5]), it is given by eq. (|6]). In terms of the problem's param- 
eters 
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Figure 8. Currents through the domain boundaries for cascade with j^n = 
Jgj as functions of time for three consecutive bursts of pair formation. The 
current flowing into the magnetosphere is shown by the soHd Hne. The cur- 
rent flowing into the NS is shown by the dashed Hne (note that the direction 
for this current is opposite to the current direction assumed in the rest of the 
paper, so it should be on average negative). Currents are normahzed to the 
Goldreich-JuHan current Jgj- The currents are averaged over 10 time steps. 



time step. These functions have spikes when secondary particles 
accelerated at the early stage of blob formation pass trough the 
corresponding boundaries. The energy carried by particles in the 
blob tail is negligible (intervals between the spikes) and so the 
energy is deposited only during bursts of pair formation. The en- 
ergy flux into the magnetosphere is larger because most of the sec- 
ondary electrons slam into the NS surface before they achieve the 
maximum possible energy while secondary positrons can gain the 
maximum energy as they fly away from the surface. The mean en- 
ergy flux going toward the NS averaged over the duration of the 
spike, e.g. over t e [6.6,8] in Fig. [9l is ~ 4 x lO^nieC^ nG]C, the 
flux going into the magnetosphere averaged over the same time is 
~ 8 X 10^ m^c^ HgjC. If the time between two successive discharges is 
(^pc/c)/ (i-e. the time between discharges is / times larger that polar 
cap flyby time) the average energy flux going into heating of the NS 
is ~ 1.5 X 10^^/"^ erg s~^cm^, this would result in the polar cap tem- 
perature Tpc ~ 4 X lO^/"^/"^ K, if the heat conductivity is neglected. 
The flux going into the magnetosphere is ~ 3 X 10^^/"^ erg s~^cm^. 
Obviously, when the time between successive bursts of pair forma- 
tion is large, the overall heating of the NS surface is significantly 
reduced. 




4 5 6 7 8 9 10 11 
t 



Figure 9. Energy fluxes trough the domain boundaries for cascade with 
jm = Jgj as functions of time for three consecutive bursts of pair forma- 
tion. The flux toward the magnetosphere is shown by the solid line; the 
flux toward the NS is shown by the dashed line. Fluxes are normalized to 
rtieC^ riojc and are averaged over 10 time steps. 

n ~ 5.6 X lOVg^^ y-j ~ 8 X 10^ , (26) 

it is still ~ 2 times larger than y^^^ given by eq. (|24] ). The par- 
ticle energy distribution at high energies is quite flat, see Fig. [6] 
So, pair producing particles are more energetic than it is expected 
from simple estimates. Because of these the total number of pairs 
generated by a single burst of pair formation should be larger than 
that assumed in "standard" Ruderman-Sutherland model. In order 
to get the final pair multiplicity detailed full cascade simulations 
are necessary, what will be subject of a separate research. The to- 
tal number of high energy particles with y > 5 x 10^ in a blob is 
~ 0.7 ^Gj^pc per cm^ of the blob perpendicular cross-section. 

In Fig. [9] I show energy fluxes trough the lower and upper 
boundaries of the calculation domain as functions of time. The 
energy flux hitting the NS surface is shown by the dashed line, 
the energy flux going into the magnetosphere - by the solid line. 
The fluxes are normalized to m^c^ Hq^c; they are computed by sum- 
ming kinetic energies of all particles leaving the domain at every 



4.5 Cascades with = 0.5, 1.5 Jgj 

When the current density in the cascade zone is diflferent from the 
GJ current density plasma flow is qualitatively similar to the case 
jm = Jgj described above. In Figs.[TOland[TT]l show snapshots of 
change density distribution and detailed characteristics of cascade 
with = 0.5 Jgj; in Figs. [12] [13]- the same plots for cascade with 
7m = 1.5 Jgj- The sizes of computation domains are the same as in 
the case of jm = Jgj. so the time t in these plots, normalized to the 
flyby time, is measured in the same units. The time intervals be- 
tween individual plots in Figs. [TO] [12] are the same as in Fig.|2]dis- 
cussed before. Note, however, that snapshots in Figs. [TT][T3l l4l are 
taken at diflferent phases of the pair formation cycle. In all cases dis- 
charges repeat quasi-periodically and creation of pair plasma goes 
through the same three phases: vacuum gap formation (? = 5.8 - 7 
for Jm = 0.5 Jgj; t = 5.483 - 6.083 for jm = 1-5jgj), formation 
and propagation of the plasma blob {t - 1 - 8.2 for jm = 0.5 Jgj; 
t = 6.083 - 7.283 for jm = 1.5jgj), and relaxation {t = 8.2 - 8.6 
for jm = 0.5jGj; t = 7.283 - 8.283 for jm = 1.5jGj). The structure 
of the plasma blob is similar - there is a charged sheath screening 
plasma in the blob from the large electric field in the vacuum gap, 
and there are large amplitude plasma oscillations in the blob. The 
superluminal plasma waves are also present. 

Cascade parameters for the case jm = 0.5 Jgj diflfer from those 
for the case jm = Jgj as follows. The size of the plasma blob is 
larger. The velocity of the previous blob's tail is smaller, Vtaii ~ 
0.44c, and the vacuum gap shrinks faster - one can see its disap- 
pearance in timeshots at ^ = 7.017 - 7.517 fFigs. [TOl[TTl) . When 
the gap closes, the charge sheath at the blob front edge disappears 
and there is no additional particle acceleration there - in Fig.[TO]at 
t - 1.461 the sheath is still present, at f = 7.567 it disappears. The 
maximum particle energy is smaller - particles do not reach the ra- 
diation reaction limited energy - but their maximum energy after 
the electric field is screened is the same yi given by eq. ( [26l ). The 
plasma density in the blob is ~ 4 times smaller. The total number 
of high energy particles with y > 5 X 10^ in the blob is ~ 0.22;igj% 
per cm^ of the blob perpendicular cross-section, what is ~ 3.5 times 
smaller than in the case jm = Jgj- The energy flux toward the NS 
is, as expected, smaller, 5.8 x 10^^/~^erg s~^cm^, and the estimated 
temperature of the polar cap is lower, Tpc ~ 3.2 x 10^/"^/"^ K. 
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Figure 10. Snapshots of charge density distribution in the calculation domain for cascade with jm = 0.5 Jgj. The same notations are used as in Fig. [21 



Cascade parameters for the case - 1.5 7gj differ from those 
for the case jm = Jgj as follows. The size of the plasma blob is 
smaller. The velocity of the previous blob's tail is Vtaii ~ 0.82c, 
it is smaller than that in the case of - Jgj* but it is signifi- 
cantly larger than Vtaii for the case jm = 0.5 Jqj- Although the gap 
shrinks faster it still leaves the domain. First generation secondary 
particles reach the radiation reaction limited energy and then slow 
down to the Lorentz factors < ji. The plasma density in the blob 
is slightly higher. The total number of high energy particles with 
7 > 5 X 10^ in the blob is ~ 0.6 Woj^pc per cm^ of the blob per- 



pendicular cross-section, what is slightly less than that value for 
jm - 7gj cascade. The energy flux toward the NS is slightly lower, 
~ 1.3 X 10^^/~^erg s~^cm^, and the estimated temperature of the 
polar cap T^, ~ 3.9 x 10^"^^"^ K. 

These diff'erences are ultimately related to the speed at which 
the tail of the previous blob leaves the domain and how many elec- 
trons are leaking from it. The cascade energetics and the ultimate 
pair multiplicity depends on the number of the first generation sec- 
ondary positrons, their maximum energy, and for how long time 
these particles are sustained at this energy. The first generation sec- 
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Figure 11. Screening of the electric field in cascade with jm = 0.5 Jgj. Snapshots are take at the same time moments as the marked snapshots in the bottom 
panel of Fig.[Tol The same quantities are plotted as in Fig. [3] 



ondary particles are accelerated up to very high energies - they can 
be accelerated up to the radiation-reaction limited energy - and are 
sustained at this energy until enough plasma is produced to screen 
the electric field in the blob. The rate of the first generation sec- 
ondary positrons production depends on how many electrons leak 
from the tail of the previous blob, and the maximum energy of these 
positrons depends on the electric field where they are injected - the 
faster the gap grow the larger the electric field. 

In the case - Jqj redistribution of particle momenta is not 
necessary for adjustment of the current density - bulk motion of 
the tail toward the magnetosphere would provide the required cur- 
rent density because the charge density is already t]qj; so in this 
case the average speed of the tail is the largest. The gap grows fast, 



and the first generation of secondary particles is created in a region 
with very strong electric field. When the current density diff'ers 
from the GJ current density, redistribution of particle momenta is 
required to sustain jm by keeping at the same time the charge den- 
sity equal to the GJ charge density. For 7m > Jgj electrons must be 
sent back to increase the current, for jm < Jgj some positrons must 
be reversed to decrease the current. The plasma as whole moves 
into the magnetosphere; low energy particles are trapped in small 
amplitude plasma oscillations and are dragged with the bulk of the 
plasma. Presence of a weak electric field would be sufficient to en- 
sure that the required number of particles on average moves toward 
the NS. Such particle reversal results in slower motion of the tail. 
When the gap upper boundary moves slower, the first generation 
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Figure 12. Snapshots of charge density distribution in the calculation domain for cascade with jm =1.5 Jgj. The same notations are used as in Fig. [21 



of secondary particles is injected in a weaker electric field, and 
the overall energetics of cascade is lower. On the other hand, the 
larger the current density the larger the number of electrons leak- 
ing from the tail; these electrons are the primary particles igniting 
the cascade. To screen the electric field at least GJ number density 
of particles is required. When the flux of the primary electrons is 
higher, the number of the first generation secondary pairs grows 
faster, and polarization of the plasma which can screen the vacuum 
electric field is achieved at smaller spacial separation; this results in 
a smaller size of the plasma blob. So, cascades with the current den- 



sity equal to the GJ current density are most energetic and should 
produce densest plasma. However, as it follows from the simula- 
tions, for jm > Jgj cascade properties seem to be less sensitive to 
the value of jm than properties of cascades with jm < joj. Hence, 
cascades with > Jq} should have energetics and final multiplici- 
ties lower but comparable to those of cascades with jm = Jgj. while 
energetics and final multiplicities of cascades with jm < Jgj should 
be significantly lower. 

In Figs.[T4]and[T5]l plot electric currents through the lower and 
upper domain boundaries for cascades with = 0.5 Jgj and jm = 
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Figure 13. Screening of the electric field in cascade with jm = 1.5 y'oj. Snapshots are take at the same time moments as the marked snapshots in the bottom 
panel of Fig. [12] The same quantities are plotted as in Fig. [3] 



1.57gj correspondingly. Except for the value of the mean current 
density these currents behave in a similar way as the currents for a 
cascade with - joj; the relative deviation of the mean over the 
cycle current density from jm is also less than ~ 10~^. 

Regarding the repetition rate of pair formation bursts in cas- 
cades with different current densities I can make only some qual- 
itative remarks. For cascades with smaller j^'s pair multiplicity is 
smaller, what must result in a less dense plasma tail; there will be 
less particles to wipe out of before the next vacuum gap can de- 
velop. On the other hand, if the current density is smaller, parti- 
cles are wiped out slower because smaller current density requires 
less particles to sustain it. For cascades with higher current den- 
sities pair multiplicity should be higher, but a larger particle flux 



is required. So, it seems that dependence of the time between the 
discharges on the current density should be moderate; it is also pos- 
sible that this dependence is non-monotonic. 



4.6 Summary of cascade properties 

The case described in details in previous sections is a good ex- 
ample of a typical Ruderman-Sutherland cascade. Insights gained 
from analyzing that case helped to draw a general physical picture. 
I performed simulation for different pulsar parameters (P, Bq, p), 
and the results are in complete agreement with the general picture 
outlined below. 

Cascades show limit cycle behavior for all physical param- 



20 A. N. Timokhin 




t 

Figure 14. Currents through the domain boundaries for cascade with jm = 
0.5 Jgj as functions of time for three consecutive bursts of pair formation. 
The currents are averaged over 10 time steps. The same notations are used 
as in Fig. [8] 
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Figure 15. Currents through the domain boundaries for cascade with jm = 
1.5 Jgj as functions of time for three consecutive bursts of pair formation. 
The currents are averaged over 10 time steps. The same notations are used 
as in Fig. [8] 



eters allowing pair creation. Pair formation is quite regular - in 
each discharge a blob of pair plasma is formed, the blob propagates 
into the magnetosphere leaving behind a tail of low energy parti- 
cles. When particle number density in the tail becomes comparable 
to Hq} a vacuum gap appears, and then a new blob of pair plasma 
forms. There are two characteristic timescales: i) time scale associ- 
ated with the size of the plasma blob Ti = Lbiob/c, ii) time between 
two successive discharges T2. The first timescale is of the order of 
Ti ~ Hrs /c for all jm, and it should be (much) smaller that the sec- 
ond timescale, ti <^T2. The simulations are inconclusive about the 
real dependence of T2 on the current density jm, but it seems that 
this dependence might be weaker that linear. 

When a new bursts of pair creation occurs, the vacuum gap 
detaches from the NS surface and propagates into the magneto- 
sphere. The upper boundary of the gap moves subrelativistically, 
while the front of the new blob moves ultra-relativistically; eventu- 
ally the gap closes. For sub-GJ current densities the tail of the previ- 
ous blob moves slower and the gap disappears faster. For cascades 
with super-GJ current densities the gap should disappears faster for 
larger values of jm, but dependence of the tail velocity Vtaii on jm 
is significantly weaker than that in cascades with sub-GJ current 



densities. The gap should persists for longest time in cascades with 

Jm — Jg] ■ 

For pulsars with large potential drop across the gap the first 
generation of pairs in cascades with jm > 0.5 Jgj reaches the ra- 
diation reaction limited energy. When electric field is screened, 
these particles propagate loosing their energy by emitting curva- 
ture photons. At first these losses are large, as radiation efficiency 
depends on particle energy as oc y^^ but then particles loose en- 
ergy more slowly. There is also a small amount of particles in the 
charge sheath at the front of the plasma blob which feel the non- 
screened electric field and are continuously accelerated. For cur- 
rent densities < 0.5 Jgj particles do not reach radiation reaction 
limited energies. However, in any case the kinetic energy of the 
first-generation particles i s larg er than it follows from estimates of 
iRuderman & Sutherland! (Il975h and so the final pair multiplicity 
and energetics of a single bust of pair formation is higher. The vast 
majority of energy is carried by the first generation pairs and so 
the heating of the NS polar cap by the cascade occurs in bursts, 
when first generation pair electrons hit the surface; heating during 
the 'relaxation' phase is negligible. If the time between successive 
discharges is large, the heating will be much lower than it is usually 
assumed in the Ruderman-Sutherland model. 

During the discharge a superluminal electrostatic wave is 
formed. As its phase velocity is larger than the speed of light it 
is not damped via Landau damping. From the performed simula- 
tions the ultimate fate of this wave is not clear because after some 
time the code stops resoling its wavelength. 



5 DISCUSSION 

I performed for the first time self-consistent time-dependent mod- 
eling of electromagnetically driven cascades which included all es- 
sential physical processes. I considered the simplest possible case 
- the Ruderman-Sutherland cascade, when all particles in the dis- 
charge zone are produced by the cascade itself. Cascade behavior is 
quite regular - spatially localized blobs of pair plasma are formed 
during regularly repeating discharges, each blob propagates into the 
magnetosphere leaving a tail of mildly relativistic plasma behind. 
Energetics of individual discharges is higher than that predicted by 
the Ruderman-Sutherland model. 

The model is one-dimensional and includes the minimum 
set of processes involved; but just because a "bare-bone" sys- 
tem was considered it was possible to get insights about gen- 
eral physics of electromagnetically driven cascades. I was inter- 
ested in dynamics of electromagnetic discharges, i.e. how particle 
are accelerated when pair production takes place. I did not fol- 
low development of the full cascade initiated by energetic par- 
ticles moving in the magnetosphere above the polar cap where 
there is no accelerating electric field. The latter problem has 
been studied before by many authors, and qualitative properties 
of cascades initiated by a relativistic particle are relatively well 
known (e.g. Daughertv & Harding 1982; Zhang & Harding 2 000l : 
iHibschman & Aronsll2001bl : lMedin & Laill2010 ). In future works I 
plan to address this problem using particle energy distributions ob- 
tained from self-consistent discharge simulations. Although quanti- 
tatively results of such full cascade simulations would be diff'erent 
from that described in the above mentioned papers, qualitatively 
they should be similar. 

In my ID simulations individual discharges are very similar 
and electrostatic oscillations are clearly visible and coherent. Usu- 
ally ID and 2D models show more coherent behavior than full 3D 
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simulations because of the inforced symmetry. Hence, the coher- 
ent behavior present in this ID simulations could be somewhat 
washed out in a more realistic 3D model. The usual picture used in 
models of polar cap discharges involves several separate discharge 
zones across the polar cap - "sparks" in terms of the Ruderman- 
Sutherland model. Whether interaction between sparks via induced 
electric field is strong enough for them to be coupled is a priori not 
clear. However, particle motion in the superstrong magnetic field in 
the pulsar polar cap is one-dimensional, the curvature of magnetic 
field lines is small, and photon trajectories only slightly deviates 
from particles trajectories; that suggests that most of the cascade 
properties deduced from current simulations can be preserved in a 
full 3D model for individual sparks, although only direct 3D simu- 
lations can prove this. 

I considered the case when particles cannot be extracted from 
the surface of the NS, but the results may be applicable to a broader 
class of problems. If cascades under diff'erent physical conditions 
work as series of discharges, their behavior should be similar to 
that described here. Namely, the blob-tail structure can be pre- 
served, pair plasma thermalization will take place, transient elec- 
trostatic wave will be excited. In particular, I suspect that polar cap 
cascades with space charge lim ited flow - in a non- stationary ver- 
sion of the model suggested bv lArons & ScharlemannI (Il979l) and 
iMuslimov & T svgan ("19921) - might have some similarities with 
the case described here. 

One of the main motivations to start this project was to study 
how cascade zone provides the current density required by the mag- 
netosphere. It has been speculated that non- stationary cascade can 
be sustained at any average current density flowing through the dis- 
charge zone (e.g. Levinson et al. 2005; Timokhin 2006). My sim- 
ulations has shown explicitly that this is indeed the case. The cur- 
rent density at any given point fluctuates strongly, but on average 
it is equal to the mean current density jm with very high accuracy. 
Current adjustment works well in the relaxation phase too, when 
both the charge and the current densities are close to the required 
values even when they are averaged over timescales much smaller 
than the flyby time of the calculation domain. It works because 
in the tail there are low energy particles trapped in small ampli- 
tude plasma oscillations, and the weak fluctuating electric field of 
plasma oscillations is able to reverse particles of both signs (at dif- 
ferent oscillation phases). Particles which flow backwards do it in 
time averaged sense, they spend more time in backward motion; 
there is no separate population of particles flowing only backwards 
all the time. Such way of current adjustment is possible because of 
efl'ective plasma thermalization which provides low energy parti- 
cles. 

The current density jm can have the sign opposite to the sign 
of the GJ charge density, jm < ( Timokhin 2006). I ran simulations 
with 7m < too; everything works exactly in the same way as in 
the case with > 0, except that the gap forms at the upper domain 
boundary, at the magnetosphere' side. This ID problem is symmet- 
ric in regard to the sign of the current density: if jm < 0, the pair 
plasma - which has positive net charge in order to sustain the GJ 
charge density - moves toward the NS and generates negative j^. 
However, at the magnetosphere' side there is no solid surface which 
can prevent charged particles from escaping; if a gap forms there, 
some charged particles can be sucked from the magnetosphere. In 
principle, it might result in presence of particles of both signs in 
the gap and, therefore, in cascade ignition at both ends of the gap; 
whether this is the case or not can not be decided based solely on 
qualitative arguments and requires accurate quantitative modeling. 
The gap will form in the tail of the blob tearing it apart; this will 



happen at distances from the NS surface larger than the polar cap 
size, and so the problem cannot be adequately described by the used 
ID approximation. The case of jm < will be addressed in a later 
work. Qualitatively, however, it seems that any cascade operating 
as a series of discharges would produce a population of low energy 
particles, and so it should be able to adjust to any current density in 
the way described above, if enough charged particles are generated. 

The simulations are inconclusive about how long the thermal- 
ization persists, because I cannot follow the plasma blob for a long 
time. It definitively works during the blob formation. For as long 
as there are low energy particles in the blob the current adjustment 
will work as described above. Note, that for current adjustment the 
number density of low energy particles should be comparable to 
Hoi, what is only a very small fraction of the plasma density in the 
blob. If, however, at some time the blob runs out of low energy 
particles, a macroscopic electric field will arise which can adjust 
the current density by creating a separate population of backward 
moving particles or/and shifting mean velocities of e l ectrons and 

iiositrons as it is suggested in e.g. IScharlemannI (Tl 91 4 ) : iLvubarskvl 
200^ . 

I performed a full-fledged kinetic modeling of pair cascades 
including all essential classes of physical processes relevant to dy- 
namics of electromagnetically driven cascades, listed at the begin- 
ning of Sec. |2] All previous attempts to model time-dependent 
cascades used on-the-spot approximation for pair injection. In 
some works fluid approximation has been used, where electrons 
and positrons were represented as fluids dLevinson et al.l I2OO5I : 
ILuo & Melrosel2008h . Although the physical situation I considered 
- the Ruderman-Sutherland cascade in the polar cap - is diff'erent 
from ones studied in previous works, it is possible to assert appli- 
cability of on-the-spot and two-fluid approximations in a general 
context, based on the general picture of cascade development in- 
ferred from my simulations. 

It turns out that the delay of pair injection due to finite time 
necessary for a photon to propagate before it is absorbed does not 
introduce new qualitative features. I also performed simulations us- 
ing on-the-spot approximation, when an electron-positron pair was 
injected at the position and at time where and when the parent par- 
ticle emitted the pair producing photon. Quantitatively, on-the-spot 
approximation introduces error in final energies of the relativistic 
particles and all depending on them cascade parameters by a fac- 
tor of several. However, qualitatively, the results are similar, i.e. 
the pattern of the plasma flow remains the same: pair formation 
is quasiperiodic with plasma blobs propagating into the magneto- 
sphere leaving tails of modestly relativistic particles. 

In the first work about modeling of time-dependent cascades 
bv lAFBeret 10(11975 1) a zero-dimensional model was used - only 
temporal, but not spatial, variations in particle number density were 
considered. In that model the production of larger amount of par- 
ticles than necessary for screening of the electric field was due 
to the time delay between the photon emissions and absorptions. 
As all later attempts to model time dependent cascades use d on- 
the spot approximations for pair injections (Levinson et al. 20051: 
iBeloborodov & Thompson! |2007I : ILuo & Melrosell2008h . the ques- 
tion about importance of pair injection delay remained unanswered. 
In my simulations the overshooting in pair number density arises 
mainly because of the spatial separation between the acceleration 
and the pair production zones in a quite regular plasma flow. Parti- 
cles are accelerated in the gap and must travel some distance before 
they can emit high energy photons. There are particles of only one 
charge sign in the gap, and so pairs are injected at only one end 
of the gap. There always exists a spatial domain with the electric 
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field (the gap) where pairs cannot be injected and the electric field 
is not regulated directly by the pair injection. The back reaction on 
the electric field proceeds only by means of gap shrinkage, which 
is slow. This causes an overshooting in pair production and so the 
intermittency of pair creation. Inclusion of spatial and temporal de- 
lays of pair injection due to photon propagation only exaggerates 
this eff'ect, but it does not introduce a new kind of behavior. Hence, 
using on-the-spot approximation in toy models seems to be justi- 
fied. On the other hand, in a situation when plasma flow can be- 
come chaotic the time delay might become a deciding factor for 
creation of plasma density overshoot. 

Two-fluid approximation, on the other hand, is inadequate. 
The pair plasma in the discharge zone acquires a large momen- 
tum dispersion and some particles become mildly relativistic. A 
weak fluctuating electric field easily reverts particles of both signs 
and plasma becomes essentially four-component (see Sec. 14.31 1. In 
two-fluid approximation at any given point at any time each parti- 
cle specie (electrons and positrons) can move only in one direc- 
tion. This i ntroduces an additiona l rigidity, which might be the 
reason why iLevinson et al ] (l2005h got strong fluctuating electric 
field throughout the whole domain. Although I considered a dif- 
ferent physical situation and the results des cribed in this paper ca n 
not be directly c ompared with the results of Levinson et al.l (I2OO5I) : 
ILuo & Melrosd (|2P08), I think that the latter are seriously flawed 
by the use of two-fluid approximation. 

Now I would like to discuss how properties of cascades could 
manifest in pulsar radioemission. Pair creation is not chaotic, with 
clear signatures of a limit cycle behavior; this ought to have strong 
observational implications. If, as it is widely accepted today, pulsar 
radioemission is directly related to pair production, the periodicity 
of cascades must be visible in power spectra of pulsar individual 
pulses. There are two characteristic time scales: ti, associated with 
the blob size, and T2, the time between discharges. The size of the 
blob to the order of magnitude is approximately the same for all 
current densities, and so ti is of the same order for any reasonable 
current density T2, on the other hand, should be more variable. 

Pulsar radioemission is highly variable on timescales compa- 
rable with the pulsar period: emission occurs mainly in form of 
subpulses, in some pulsars subpulses drift, some pulsars changes 
modes and/or switches ofl" for many periods. This hints that cur- 
rent density can fluctuate because of some processes in volving the 
whole magnetosphere (e.g.'Arons*1983l: lTimokhinll201Q ). Cascades 
can adjust to any reasonable current density, and so the current den- 
sity at a fixed colatitude might vary on timescales much larger that 
Ti , T2; on the other hand, the current density varies across the pulsar 
polar cap anyway. Because of these, an individual subpulse repre- 
sents emission averaged over time and space, or, in other words, 
over a range of diff'erent jm's, and so the features of cascades along 
separate field lines will be smeared. Hoverer, the time scale asso- 
ciated with the size of plasma blob ti by the order of magnitude 
remains the same and should be clearly visible in the power spec- 
trum. The second time scale T2 should be less prominent, but, as 
discussed before, it might be not extremely sensitive to the current 
density, and, therefore, it might manifest as a broad feature in the 
power spectrum. 

The blob is of the same length as the region with the accelerat- 
ing electric field. In the Ruderman-Sutherland model this length is 
small, and the corresponding timescale Ti is less than a microsec- 
ond. In space charge limited flow models, on the other hands, the 
length of the acceleration zone should be comparable to the NS ra- 
dius; if in this case cascades work similarly, ti should be of the 
order of ~ 100 //sec. The second time scale, T2, should be sub- 



stantially longer, a factor from few to hundreds. There are evi- 
dences of difl'erent characteristic time scales in pulsar microstruc- 
ture, from nano- to milliseconds; in some pulsars microst ructure is 
also q uasiperiodically modulated (e.g Boriakoff 1976; P opov et al.l 
I2OO2I) . It is not clear whether microstructure timescales are due to 
polar cap cascades variability or not, but ti , T2 can be in the range of 
observed microstructure modulation times, and cascades operating 
as a series of discharges should have double-timescale signature. 

The problem of pulsar radioemission mechanism in notori- 
ously difficult and currently there is no reliable theory which could 
adequately explain it. The firmly established observational fact 
about pulsar radioemission is that it is due to some collective pro- 
cess. In my simulations I saw formation of a large amplitude elec- 
trostatic wave. Its phase velocity is larger that the speed of light, 
and it is not damped via Landau damping. In one dimension in a 
superstrong magnetic field only electrostatic waves exist, but in a 
real pulsar such wave can be coupled to an electromagnetic mode; 
if it stays superluminal, it can escape the magnetosphere. This wave 
is a collective form of emission, as it involves coherent macroscopic 
plasma motion. The simulations are inconclusive about the fate of 
that wave because at some point the numerical scheme stops re- 
solving its wavelength; it is also not clear how coherent the whole 
picture is in 3D. May be it is too preliminary to tell whether pulsar 
radioemission, or some of its component, is related to this transient 
wave, but in future research special attention should be paid to such 
transient waves. 
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APPENDIX A: ONE-DIMENSIONAL TIME-DEPENDENT 
ELECTRODYNAMICS OF THE POLAR CAP 

In t he reference frame corotating with the NS the Gauss law is (see 



in t he reterence trame corotating 
e.g. lArons & Scharlemamill 19791) 

V X £ = 47T{7] - T/Gj) . 



(Al) 



In the ID approximation the only changing component of electro- 
magnetic fields is the electric field parallel to the static magnetic 
field of the NS. The solution of equation dAll) is given by 
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Figure Al. Current flow through the surface on the NS. See text for expla- 
nation. 



E = El=o +47T I (77 - t/gj) ds 
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(A2) 



I am solving a non- stationary problem where boundary conditions 
can change with time because the magnetosphere can response to 
the changes of conditions in the polar cap. If the electric field just 
outside the NS surface E\^=q is known at any given moment of time, 
the electric field in the calculation domain can be calculated using 
eq. (IA21) . For this problem it is more convenient to reformulate the 
boundary conditions on the electric field at the NS surface E\^=q in 
terms of the electric current flowing through the system. 

As 7/gj does not change with time, diff'erentiating eq. lA2l with 
respect of time and using charge conservation 



dt dx 
I get 



dE^ 

'dt 



dE 
'dt 

or 



where 
, _ I dE 
4^^ ^ 
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(A5) 



(A6) 



To clarify the meaning of let us consider a small region at the 
NS surface, see Fig. lAll NS crust can be considered as a good 
conductor; the charges can accumulate only on its surface, and the 
electric field in the crust is zero. The electric field at the NS surface 
E\^=Q = 47TO-, where <t is the surface charge density. The change 
of the total charge in the fiducial volume in Fig. lAll Sq is due to 
currents through the boundaries of the volume: 



,6S 

For the electric field at the NS surface I have then 
1 dE 
47r 'dt 

Substituting this expression into eq. ( IA6I ) I get 



dcr dJin 



(A7) 



(A8) 
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;™ = ^, (A9) 

dS 

i.e. jm is the current density which flows in the NS crust toward 
the discharge zone; it causes current in the discharge zone and/or 
accumulation of charges at the NS surface. In other words, is 
the current density that the magnetosphere wants to flow through 
the cascade zone. Eq. ( IA5I ) is a convenient form for an equation 
for the electric field in a problem where a large system with very 
high inductivity requires some specific current density from a much 
smaller systeni plugged into the same electrical circu it (see e.g. 



Levinson et al.l l2005inBeloborodov & ThompsonI 12007 1). Note that 



eq.|A5] correctly accounts for the retardation of changes in the elec- 
tric field - at any given point in space the electric field changes if 
j deviates from the current density j is generated by particle 
motion, and the latter cannot move faster than the speed of light. 



